#pragma once
#include "gauss.h"
#include <math.h>

// 定义一个需要积分的函数
double myFunc(const std::vector<double> &x)
{
    double y = 0.0;

    y = 3.0 * pow(x[0], 2.0) + x[1] + pow(x[2], 4.0);

    return y;
}

double gauss::getIntegration()
{
    std::vector<double> x;

    double sum = 0.0;
    for (int i = 0; i < gaussPointNum; i++)
    {
        for (int j = 0; j < gaussPointNum; j++)
        {

            for (int k = 0; k < gaussPointNum; k++)
            {
                x = {kesi[i], kesi[j], kesi[k]};
                sum += w[i] * w[j] * w[k] * myFunc(x);
            }
        }
    }

    return sum;
}

gauss::~gauss()
{
}

gauss::gauss(const int &i)
{
    gaussPointNum = i;
    kesi.resize(i);
    w.resize(i);
    switch (i)
    {
    case 1:
        kesi[0] = 0;
        w[0] = 2;
        break;
    case 2:
        kesi[0] = -1.0 / sqrt(3.0);
        kesi[1] = 1.0 / sqrt(3.0);
        w[0] = 1;
        w[1] = 1;
        break;
    case 3:
        kesi[0] = 0;
        kesi[1] = sqrt(3.0 / 5.0);
        kesi[2] = -sqrt(3.0 / 5.0);

        w[0] = 8.0 / 9.0;
        w[1] = 5.0 / 9.0;
        w[2] = 5.0 / 9.0;
        break;
    case 4:
        kesi[0] = sqrt((3 - 2 * sqrt(6.0 / 5.0)) / 7.0);
        kesi[1] = -sqrt((3 - 2 * sqrt(6.0 / 5.0)) / 7.0);
        kesi[2] = sqrt((3 + 2 * sqrt(6.0 / 5.0)) / 7.0);
        kesi[3] = -sqrt((3 + 2 * sqrt(6.0 / 5.0)) / 7.0);

        w[0] = 0.5 + 1.0 / (6.0 * sqrt(6.0 / 5.0));
        w[1] = 0.5 + 1.0 / (6.0 * sqrt(6.0 / 5.0));
        w[2] = 0.5 - 1.0 / (6.0 * sqrt(6.0 / 5.0));
        w[3] = 0.5 - 1.0 / (6.0 * sqrt(6.0 / 5.0));

        break;

    default:
        std::cout << "Guass Integration Point Num not in Range [1-4]:" << std::endl;
        std::cout << "Press Enter to Stop!!" << std::endl;
        system("pause");
        exit(0);
        break;
    }
}
